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ABSTRACT 

We present a new implicit numerical algorithm for the calculation of the time depen- 
dent non-Local Thermodynamic Equilibrium of a gas in an external radiation field that 
is accurate, fast and unconditionally stable for all spatial and temporal increments. 
The method is presented as a backward difference scheme in 1-D but can be readily 
generalised to 3-D. We apply the method for calculating the evolution of ionisation 
domains in a hydrogen plasma with plane-parallel Gaussian density enhancements 
illuminated by sources of UV radiation. We calculate the speed of propagation of ion- 
ising fronts through different ambient densities and the interaction of such ionising 
fronts with density enhancements. We show that for a typical UV source that may be 
present in the early universe, the introduction of a density enhancement of a factor 
~ 10 above an ambient density 10~^ cm~'^ could delay the outward propagation of an 
ionisation front by millions of years. Our calculations show that within the lifetime of 
a single source a few million years), and for ambient intergalactic densities appro- 
priate to redshifts z ^ 6 — 20, degrees of ionisation of ~ 10"'^ — 10~^ can be achieved 
within its zone of influence. We also present calculations which demonstrate that once 
started, ionisation will proceed very efficiently as multiple sources are subsequently 
introduced, even if the time between the appearence of such sources may be much 
longer than their lifetimes. 

Key words: numerical radiative transfer, early universe, proto-galactic clouds. 



1 INTRODUCTION 

In local thermodynamic equilibrium it is explicitly assumed 
that the local kinetic temperature of the gas determines the 
level populations of the atoms or molecules via the Boltz- 
mann distribution. This situation pertains when collisions 
dominate or when the radiation field is a black body cor- 
responding to the local kinetic temperature of the gas. In 
contrast, in non-LTE, the level populations are determined 
by the rate equations, and depend both on the local kinetic 
temperature and on the local radiation field which may not 
be black body like (see Oxenius (1986) for precise defini- 
tions). In most problems in astrophysics, non-LTE is asso- 
ciated with steady state level populations. 

There are situations where time dependent non-LTE 
calculations are required. This could occur for instance when 
the hydrodynamical time scales are shorter or comparable 
to the radiative time scales (e.g. in the formation of shock 
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fronts or in stellar explosions). Another situation of current 
interest where time dependent effects are important occurs 
in the description of the re-ionisation of the universe (e.g. 
Peacock (2000), Liddle and Lyth (2000)). Here, the densities 
are sufficiently low and the recombination times correspond- 
ingly large that a global equilibrium between photoionisa- 
tion and recombination may not be achieved during the life- 
time of the source of the ionising photons. The modelling 
of this phenomenon in an inhomogeneous medium with a 
complex 3D geometry is a formidable computational task. 

The past few years have seen the development of many 
numerical methods aimed at addressing the problem of re- 
ionisation in a 3-D framework. Rather than attempting to 
solve this problem in its full complexity, many of these meth- 
ods make practical approximations to the solution of the 3-D 
transfer equation. We note in particular adaptive ray trac- 
ing methods (Abel, Norman & Madau 1999, Abel & Wan- 
delt 2002, Sokasian et al. 2003), the fast Fourier transform 
method (Cen 2002) and the moment equation method us- 
ing Eddington tensors with an optically thin approximation 
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(Gnedin &Abel 2001). In all of these methods, a key ele- 
ment is the coupling between the chemical (rate) equations 
and the radiative transfer equations which determines how 
well the ionisation front is resolved. In general no explicit 
allowance is made for the time rate of change of the spe- 
cific intensity (the term) in the transfer equation so 
that the methods yield faster than light propagation of ion- 
isation fronts close to the source. Although these methods 
cannot be rigorously justified in all their assumptions, they 
are physically motivated, and provide a description of the 
process of re-ionisation at some level. Another approach has 
been to use Monte Carlo methods (e.g. Masselli, Ferrara & 
Ciadi 2003, Masselli et al. 2004). The major limitation here 
is that they require the adoption of an explicit time step- 
ping method which is unstable for large time steps. This, 
combined with the large number of trials that are needed, 
appears to limit the method by the computational efi'ort that 
is required. 

In this paper we present a new and robust numerical 
method for solving the time dependent rate equations and 
the time dependent radiative transfer equation concurrently 
using a backward differencing scheme. In section 2 we dis- 
cuss the nature of the problem and present our numerical 
scheme formulated in 1-spatial dimension. Illustrative re- 
sults describing the time evolution of ionisation domains 
that propagate through homogeneous and inhomogeneous 
plane parallel distributions of initially neutral hydrogen gas 
are presented in section 3. The effects of including multiple 
sources are discussed in section 4. The models in these two 
sections can serve as a useful bench mark for comparisons of 
results from other codes and for understanding the process 
of re-ionisation in the early universe. Our main results are 
summarised in section 5. 



2 THE BASIC EQUATIONS AND 
NUMERICAL METHOD 

We consider only the simplest case of the re-ionisation of 
hydrogen although the method can be easily extended to 
include helium. The basic equations are then the rate equa- 
tions which describe the populations of neutral hydrogen 
(HI) and ionised hydrogen (HII), and the radiative transfer 
equation which describes the evolution of the mean radi- 
ation intensity due to the absorption and the emission of 
radiation. 



2.1 Formulation for a two level atom 

In the simplest case that we consider, the model H atoms 
have one bound and one continuum state. The two processes 
of relevance are then ionization from the lower (ground) 
state and recombination from the upper (continuum) state. 

The rate equation for the number density of neutral 
hydrogen atoms Ng can be written as 

dNoix,t) 



dt 



-aJ{x,t)No{x,t) + a{Ntot - Noix,t)y (1) 



where A^'tot = No + Ni and A^i is the number density of 
ionised hydrogen (HII). Here a is the photoionisation cross 
section, a is the recombination coefficient, and 
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Figure 1. Sequence of steps to be used to cover the x,t plane. 
The starting position is at the lower left corner. The thick lines 
and large dots indicate the where the initial condition are to be 
specified. 




Figure 2. The convergence to the stationary state. The curves 
show from top to bottom the evolution of the Model b3 (see 
Table 1) density distribution after 0, 2, 4, 8 Ma (assuming that 
the illuminating source lives for this period) and the correspond- 
ing steady state distribution calculated by a completely differ- 
ent algorithm (see text). For clarity the latter is shifted by 
AlogTVo = —0.5, otherwise the two lowest curves would overlap. 



is the mean intensity of the radiation field where du) is the 
element of solid angle. Note that all occupation numbers are 
functions of position x and time t. Since the radiative rates 
depend on the local mean intensity, it is necessary also to 
solve concurrently the radiative transfer equation which we 
give here for a ray in the direction n — (1, 0, 0) 



ldl{x,t) dl{x,t) 



dt 



+ 



da 



= ^aNo(x,t) {I{x,t) - S{x,t)l3) 



J{x,t) 



J(x, t, n)dLj 



(2) 



where S{x,t) is an appropriate source function. Equations 
and 121 have to be solved as an initial value problem with 
7Vo(a;,0) — Ntot{x), Ni{x,0) = with the inflow boundary 
condition l{0,t) — Urr on the radiation field where Urr is 
specified. That is, we assume that the medium is completely 
neutral at the time t = Q when the source of ionising photons 
first turns on. 

To our knowledge there does not exist an analytical so- 
lution for equation combined with equation |3 which is 
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Table 1. Basic model data 



Model 


^ *amb 


*enn 








"^start 


^1 1^3 


At 




_ "3 

cm 


cm-3 


10^4 cm 


10^1 cm 


lO'' erg cm— ^ s— 4 


104** cm-2 


10^ cm s— 1 


1043 s 


a 


10"^ 









1.0 


0.00064 


1.3 




b 


10-4 









1.0 


0.0064 


0.96 




c 


10-3 









1.0 


0.064 


0.20 




d 


10-3 









3.0 


0.064 


0.30 




bi 


8.61 10-5 


9.47 10-4 


31 


0.5 


1.0 


0.0064 


1.1;1.2;0 


8 




7.83 10-5 


8.61 10-4 


31 


2.0 


1.0 


0.0064 


2.0;1.7;0.07 


26 


bs 


4.74 10-5 


4.93 10-4 


31 


4.0 


1.0 


0.0064 


1.8,2.5;0.13 


31 


b4 


4.94 10-5 


5.12 10-4 





4.0 


1.0 


0.0064 


0.55;0.67 




si(two sources) 


10-4 


10-3 




4.0 


1.0 


0.0064 







a system of non-linear partial differential equations. There- 
fore, we have to solve the discretized equations numerically. 
An immediate possibility is to employ the method of lines 
(cf. Schiesser, 1991) which would incorporate an up-wind 
discretisation in the spatial coordinate and the use of a solver 
for the resulting system of ordinary differential equations (as 
e.g. a Stoer-Burlisch method). Unfortunately, since the right 
hand sides of both the rate and the transfer equations are 
always negative in our cases, such an discretisation leads to 
instabilities whenever the temporal stepsize exceeds a cer- 
tain value (cf. Kahaner et al., 1989). Moreover, the system 
is often very stiff as a consequence of the large differences in 
the time-scales associated with photoionisation and recombi- 
nation. As a consequence, in explicit schemes unacceptably 
short time steps have to be employed. Therefore, we em- 
ploy here an implicit scheme, i.e. an up-wind discretisation 
in x-t space. Only first order differencing is considered here 
in order to avoid overshooting. In general, this involves the 
solution of a huge system of simultaneous non-linear equa- 
tions (cf. Schiesser, 1991). Fortunately, this is not necessary 
here, when the fact is exploited that light rays propagate 
always in one direction, i.e. the specific intensity at a spa- 
tial point x^^^ depends only on the specific intensity at 
and the absorption coefficients aNo and the source function 
S m dx = x^ . . .x^'^^. If we restrict ourselves to linear ap- 
proximations these quantities can be assumed to be given 
by aNo{x^+^,t) and S'(a;'+\t). 

The discretized equations and |21 therefore read 

A^o(a:'+\t''+i) - Afo(a:'+4,tfc) 

= ~(jNo(x'+\t''+^)I{x'+\t''+'') 

+a{N,,,{x'+') - No{x'+\t''+^)f (4) 

I{x'+^,t''+^)- I{x\t''+^) 
1 J(3:'+\t''+i)-J(a:'+4,ffe) 

= -GNo{x'+\t''+^) 

X \l{x'+\t''+^) - S(x*+\ *''■+')] (5) 

or if we solve equation |3 for /(x'"'"'^, t*"'"^) and insert the 



result into Eq.0]we get 

Aro(a:'+4,t''+4) - Afo(a:'+4,tfc) 

= -crAfo(2;"+4,t'=+i) 

xi \l{x'+\t''){x'+^ - x') + I{x\t''+^)c{t''+^ - t'') 
+aNo{x\t''+^)c{t''+^ - t''){x'+^ - x')Six'+\t''+^)^ 
+a (7Vt„,(x"+4)-iVo(x'+\t'=+^))' 

(6) 

and 

I{x'+\t''+^) 

= i [/(x'+\t'=)(x'+' -X*) +/(x\t'=+i)c(t'=+i -t'') 
+aNoix\t''+'')c{t''+^ -t''){x'+'' ~x')S{x'+\t''+'')^ 

(7) 

with 

d ^ {x'+'' - x') + c{t''+^ - t'') 

+ (x'+'^ -x')c{t''+'' -t'')aNo{x'+\,t''+^) (8) 

Equation 1^ is a cubic equation for No{x'^'^^ ,t'''^^), i.e. it is 
an equation in one variable, whenever I{x^,t^'^^) is known. 
Therefore, if we make use of the boundary conditions, we can 
march through the x — t plane as indicated in Figure The 
equation can be solved analytically or by a simple numerical 
solver. This is much more economic in memory and CPU 
requirements than the general case, where one has to solve at 
every time-step a non-linear system, in which each equation 
depends on the densities at all spatial gridpoints. In this 
way we get an algorithm that is unconditionally stable for 
all spatial and temporal stepsizes and that is very efficient 
even for very many grid points. 

We demonstrate that our method converges to the re- 
quired steady states for Model b3 (see Table 1) in Fig.|21 In 
order to calculate the steady state model by an independent 
method we have first set dNo/dt = in the rate equation, 
solved the algebraic equation for A'^o and inserted the result 
into the transfer equation which we subsequently solved as 
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Figure 3. The time evolution of the logarithm of the num- 
ber density of neutral hydrogen No{x,t) at three different times 
t = 0, 2, 4 Ma in uniform density models. From top to bottom: 



tot 



10" 



cm 
= 10^ 



erg cm" 



^ erg cm~2 (Model a); Ni 



cm 

10^ erg cm-2 (Model c); A^tot 
cm -2 (Model d) 



10 

2 (Model b); JV, 
10-3 



tot - 
3 



= 10" 



tot = 10"^ 

cm" 3, 

= 3 X 10^ ei 



an ordinary differential equation by means of a predictor- 
corrector scheme. 

Depending on the computer and operating system as 
well as the software used the CPU times range from a few 
seconds to a few minutes. This is about a factor 1000 more 
efficient than an explicit method of lines approach. 

2.2 The inclusion of many atomic levels 

If the details of the ionising flux are to be taken into account, 
the relevant wavelengths A and the wavelength dependence 
of opacities and of the source function have to be consid- 
ered explicitly. Therefore the rate equation and the transfer 
equation read 



dNo{x,t) 
dt 



-No{x,t) / a{X)I{x,t,\)dX 
Jo 

n2 



+a{Ntot ~ No{x,t)y 



(9) 



1 dI{x,t,X) ^ dI{x,t,X) 
c dt dx 

= -a{X)No{x,i) {I{x,t,X) - S{x,t,X)) 

(10) 

The discretized versions (cf. equations |^ can now be 
written 

No{x'+\t''+'') - No{x'+\t'') 

L 

-L.(^I(x^+\t\x'){x'+'^x^) 
+I{x\t''+\x')c{t''+^ 

+a{X')No{x\t''+'-)c{t''+^ - t''){x*+^ - X*) 
xS{x'+\t''+\x') 



+a (^Ntotix'+^) - Noix'+\t''+')y 



(11) 



with 



d(X') = {x'+^ - x') + c{t''+^ 



r)o-(A')Afo(a;"+\r+') (12) 



(w^ are the weights for the wavelength integration). It is seen 
that the first equation contains only No{x'^^ ,t'''^^) as an 
unknown. However, the polynomial for the determination of 
No{x'^^^ ,t'''^^) is now of degree 2L+1. Fortunately, this does 
not slow down the calculations significantly if a numerical 
solver is used. 



2.3 Generalisation to more than one source of 
radiation 

In the general case, we need to consider the presence of more 
than one source of radiation. We describe the algorithm we 
use in terms of two sources of radiation, but the generalisa- 
tion to many sources is immediate. 
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Suppose we have two sources of illumination: one is sit- 
uated at s = and shines from time t = to time t — ti 
with constant luminosity and the other one is placed at 
X = Xend and shines from t = t2 to t = t^nd where the 
subscript 'end' indicates the maximum value. We now em- 
ploy the two-stream-approximation where Ip indicates the 
direction of increasing x values and Im the opposite direc- 
tion. The mean intensity is then given by J = {Ip + Im)/2. 
Corresponding neutral hydrogen densities are denoted by 
Nop and Nom- The solution is derived now by a fixed-point 
iteration for No{x,t). For the first step l}n{x,t) = is as- 
sumed for all x and t values. Ip{x,t) and NQp{x,t) are then 
calculated for all spatial and temporal grid points as de- 
scribed above. Next, I^{x,t) and A'()„(a;,f) are calculated 
starting from t — and x = Xend (using the same scheme 
as above but just going in the opposite direction) with the 
Ip{x,t) as just derived. NQ{x,t) is then calculated as the al- 
gebraic mean of NQp{x,t) and NQ^{x,t). Next, Ip(x,t) and 
N§p{x,t) are evaluated with the assumption that the inten- 
sity in the opposite direction is given by the values of the 
previous iteration, and then and NQ^{x,t) leading to 
NQ{x,t). The next iteration then follows, and we stop the 
process when \NS'{x,t) - N^''^ {x,t)\/NS' [x,t) < e for all 
spatial and temporal grid points. For e = 10~^ this was al- 
ways the case for m ^ 5 in our computations, i.e. we find a 
fast convergence. 

For the models discussed in this paper, the term 
is unimportant and is therefore not included. The role played 
by this term is discussed in section 3.3. Since we are not 
dealing with recombination radiation in this paper, we set 
S{x,t) = in our calculations. Furthermore, unless explic- 
itly stated, all our illustrative models assume two levels. 



3 MODEL RESULTS FOR SINGLE SOURCES 

3.1 Plane parallel uniform density slabs 

We assume that the medium is plane parallel extending from 
X = to X = a^max. The density profile is specified to be 



iVtot(a;, 0) = iVamb + iVonh exp[-(- 



In 



(13) 



where A^amb is the ambient density. The second term allows 
for a density enhancement with a Gaussian profile of ampli- 
tude A'cnh and half width 1^ centered at x = Xa. 

We choose values for the free parameters to correspond 
approximately to the conditions that are expected in the 
universe when the first sources of radiation turned on. For 
the standard cosmological parameters, the gas density at 
redshift z is given by iVtot ~ 10-^(1 + zf so we consider a 



density range 10 



10 cm appropriate to z 



6-20. 



We set the density enhancements to peak at ~ 1 dex above 
the mean. 

We place the source of UV photons at a; = and ob- 
tain results for a spatial region that extends to Xmax ~ 20 
kpc. The inflow intensity at the outer boundary of the plane 
parallel region that is being irradiated will depend on the lu- 
minosity of the source, its geometry, and its location relative 
to the region. In our examples we use an incident intensity 
Jir,, = 10^ erg cm~^ s~^ which corresponds to the total in- 
tensity of a black body of temperature T — 10^ K (a UV 



source) diluted by a factor ~ 10"^". Since these are 1-D 
models, we do not specify the nature of the sources or of the 
dilution. 

Equation can be used to estimate a recombination 
time scale 

No 



^{Nu 



No 



(14) 



Initially, when the gas is mainly neutral, recombination is 
ineffective {ticc oo), and the time evolution is governed 
by photoionisation. As the degree of ionisation increases, 
tree decreases, and for a gas that is 50% ionised, the above 
expression yields a characteristic recombination time scale 
(using a ~ 10"^** cm^ (Osterbrock (1974)) 



6.3 X lO'' 



rlO~ 



iVt< 



Ma 



(15) 



This time scale becomes smaller, the higher the density. 

Likewise, using a ^ 10~^* cm^ (Osterbrock 1974) equa- 
tion |H] yields a photo- ionisation time scale 

1 ,10^ erg cm^^s"-^, , , , , 

tphot 7 ~ 0.32 [ ^- ] Ma (16) 

This time scale is smallest for the highest intensities. As the 
intensity of radiation decreases (due to photo-absorption), 
the photo ionisation time scale increases, with tphot ^ oo 
as / ^ 0. In principle, the radiation held and the ionised 
fraction can adjust so as to reach equilibrium at all points 
in the medium. However, in the situations that we consider 
the level populations never reach full equilibrium during the 
lifetime of the source. Our calculations, however, allow us 
to chart the rate at which this equilibrium is approached at 
different times at different points in the medium. 

In the first set of calculations we do not allow for a 
density enhancement (A'cnh ~ 0) and consider three models 
with different uniform densities A^amb = 10~^ cm~^ (Model 



A^a, 



10--* cm-^ (Model b) and A^amb = 10" 



(Model c). Figure 121 shows the run of density of neutral hy- 
drogen for these three cases at three different times t = 0,2,4 
Ma. For the model with the lowest ambient density (Model 
a), the entire computational region is optically thin to UV 
radiation (even a.t t — 0) so that the effects of radiative 
transfer are minimal. As a consequence, the density of neu- 
tral gas is nearly independent of x, but decreases with 
time. The decline occurs exponentially on the photoionisa- 
tion time scale being dominated by the first term in the 
rate equation Q The neutral component does not show any 
marked spatial structure at a given time within the domain 
of calculation up to the time tmax when the calculation is 
terminated. 

As the ambient density is increased by a factor 10 
(Model b; Figure |3J optical depth effects become more im- 
portant. Regions that are closest to the source are least af- 
fected by transfer effects and ionise more quickly in com- 
parison to regions that are further away. As a consequence 
recombinations become relatively more important as one 
moves away from the source. This is shown by the grad- 
ual increase in the neutral gas density with increase in x at 
any given time t > 0. In this model there are clear changes 
in which - as we shall see - can be used to define the 

dt 

location of the ionisation front. 

As the ambient density is increased even further (Model 
c; Figure|2l, optical depth effects dominate except very close 
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L |Hl| 



1 llc^t 




Figure 4. The rate of change of the neutral density l^^^pl on a logarithmic scale for Models a (top left),b (top right), c (bottom left) 
and d (bottom right). The front corresponds to the region of highest gradient shown here with the lightest shading. The speed of the 
front is seen to decrease as the density increases for a given incident intensity, and to increase as the incident intensity increases at a 
given density. 



to the source. A well defined ionisation front which prop- 
agates outwards with time is now clearly seen within our 
computational domain. The effect of increasing the incident 
intensity by a factor of three to 1^^,. = 3 x 10^ erg cm~^ s"'^ 
for the same density distribution is shown in Model d. The 
ionsiation front is now seen to propagate further in the same 
time interval. 

Comparing the models in Figure^ we see that the max- 
imum degree of ionisation that is achieved at t = tmax ranges 
from 10"^ to 10"^ which agrees in order of magnitude with 
the mean observed degree of ionisation of the difluse inter- 
galactic medium at the tail end of re-ionisation (Fan et al. 
2003). 

The properties of the fronts are best illustrated in Fig- 
ure ^ where we present | \ as a function of x and t for 
Models a, b ,c and d. We identify the front as the region of 
the maximum rate of change of ionisation (in an absolute 
sense). These diagrams clearly show how the front propa- 
gates in time and changes its spatial structure. The proper- 
ties of the fronts, such as the speed, can be extracted from 
this data and these are given in Table 1. 

In Model a, the front propagates fast as it eats its way 



through the low density medium. The speed of the front 
is highest for this model. In Model b, the higher density 
slows down the front, so that it advances at a slower rate 
while in Model c, the even higher density reduces the speed 
further. We note that in these calculations, the front speed 
is essentially independent of time except close to t = where 
the speed is not well defined (see section 3.3). 

The front speed is a strong function of the irradiating 
intensity as shown in Figure 2] In general, a stronger source 
will result in a front that propagates faster qualitatively sim- 
ilar to a model with a lower mean density. 

The spatial distribution and time evolution of the inten- 
sity of radiation for Models a,b and c are shown in Figure 
|S] As already noted, this sequence corresponds to one of in- 
creasing optical depth at time t = when the gas is neutral. 
Very close to the source, the intensity is maintained at essen- 
tially the input intensity at all times considered due to the 
dominance of photoionisation. The intensity that emerges 
from the computational region (i.e. at a; = a;max), on the 
other hand, evolves with time at different rates depending 
on the density of the model. In Figure |K| (Model a) , the 
emergent intensity is essentially the input intensity with lit- 
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log I 
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S t lHa.] 



Figure 5. The spatial distribution and time evolution of the 
logarithm of the intensity of radiation of Models a (top left), b 
(top right) and c. These calculations show how the source emerges 
from the dark phase to become bright and visible as the front 
propagates through the computational volume 



tie change as t approaches t = tmax since photoionisation 
dominates over the entire computational volume. On the 
other hand, in Figure |3 (Model b), this intensity increases 
from a low value at initial times to a value close to the 
input intensity as t approaches imax and a larger fraction 
of the gas becomes ionised and therefore less opaque. The 
source will therefore first appear dimmed when viewed from 
^ = a^max brightening as t increases towards fmax- In a sense, 
this source emerges from a dark phase (the 'dark ages' in a 



-6 
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X [kpc] 




10 

X [kpc] 




20 



X [ kp c ] 



5 10 15 20 

X [kpc] 

Figure 6. The time evolution of the logarithm of the number 
density of neutral hydrogen No{x,t) at three different times t = 

0, 2, 4 Ma in a non uniform medium with density enhancements 
described by Gaussian profiles. The parameters of the four models 
(Models bl, b2, b3, b4 from top to bottom) are given in Table 

1. All models are normalised to have a total column density S = 
6.4 X 10^* cm~^ corresponding to Model b 
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cosmological context) within the maximum time used in our 
computations. In contrast, at the even higher densities rep- 
resented in the model in Figure |3 (Model c), the bulk of the 
material remains neutral for all t ^ tmax so that the source 
never emerges from the dark phase. Clearly, had these calcu- 
lations been carried out for even larger times, the ionisation 
front will eventually move outwards to a;max and the source 
will no longer be dark. We note that apart from a scaling 
factor, the inverse of the intensity in Figure|5]is proportional 
to the photoionisation time scale. 

3.2 Plane parallel slabs with Gaussian density 
enhancements 

In the second set of models we introduce a density enhance- 
ment in the middle of the computational region at Xs — 10 
kpc (Models bl, b2 and b3) and at = (Model b4) 
parametrised by N^rih and The widths of the enhance- 
ments increase from Model bl through to Model b3, and 
the densities are normalised so that the total column num- 
ber density is always E = 6.4 x 10^* cm~^. The parameters 
for this series of models are also summarised in Table 1. 

We show in Figure|S|a, b, c and d the density of neutral 
hydrogen for Model bl, b2, b3 and b4 at three different times 
f = 0,2,4 Ma. 

The introduction of a density enhancement leads to a 
local increase in the extinction coefflcient and therefore a 
decrease in the number of ionising photons that can propa- 
gate to the far side of the enhancement in Models a, b and c. 
As a result ionisation is more effective on the 'near side' of 
the density enhancement (the side that faces the source of 
irradiation) than on the 'far side' as can be seen by the drop 
in the level of ionisation as one crosses the enhancement. 
The model with the largest value of 1^ (Model bS) shows 
the strongest drop in the level of ionisation at any given 
time. As t increases, the overall level of ionisation increases 
and the optical depth through the density enhancement de- 
creases. The relative change in the level of ionisation across 
the enhancement also decreases. 

In Model b4, the source is at the peak of the density 
enhancement so that the radiation propagates through a 
density distribution that decreases with distance from the 
source. The time evolution is very similar to what is found 
for the far side of the density enhancements in Models a, 
b and c. The ionisation proceeds slowly first close to the 
source, and then propagates outwards at later times. 

The interaction of the outward propagating front with 
a density enhancement is beautifully illustrated in Figure 
[7| which is analogous to Figure 2] Here we see that as a 
front with a velocity vi encounters a density enhancement, 
it slows down dramatically to a mean speed V2 as it adjusts 
to the higher density environment. Eventually the front pen- 
etrates the density enhancement and emerges at a speed 113. 
Ill and fs correspond very approximately to the speeds ap- 
propriate to the ambient density and intensity on the near 
and far sides of the density enhancement. However, in gen- 
eral V3 is larger than vi with the difference being largest for 
the enhancement with the largest Gaussian half width. The 
encounter creates a well defined delay At/ in the process 
of ionisation which is a function of the parameters of the 
Gaussian. 

In the case of Model b4, where the source is at the peak 



of the density enhancement, the velocity is initially low, and 
then steepens very quickly as the front breaks through and 
the velocity approaches the value appropriate to the ambient 
density. 

The time delays and the other properties of the fronts 
are summarised in Table 1. 

The above calculations show that in the presence of 
Gaussian density enhancements, even of modest amplitude, 
can result in significant delays in the onset of re-ionisation. 
For suitable choices of parameters (as in Model bS) that may 
be appropriate to the first sources of radiation that form in 
the early universe, the delay could exceed the lifetime of the 
source. Once the source dies out, the gas would recombine 
on the much larger recombination time scale free. Pockets 
of ionised regions would thus survive enshadowed by high 
density regions until another source of UV photons turns on 
in its vicinity and proceedes to ionise its surroundings. 

The spatial distribution and time evolution of the in- 
tensity of radiation for the Models bl, b2 and bS are shown 
m Figure H These fi gures should be compared with the cor- 
responding uniform density Model b in Figure |S] which has 
the same total column number density. The intensity is seen 
to drop sharply across the density enhancement initially due 
to the local increase in optical depth, but with time, as the 
region ionises and becomes less opaque, the intensity ap- 
proaches the input intensity. 

We conclude this section by presenting calculations of 
the recombination time scale and the actual time scale on 
which the number of neutral hydrogen atoms change in 
Model b3 (Figure |HJ as a function of {x,t). The diagram 
shows that initially the recombination time scale has a min- 
imum at the peak of the Gaussian density enhancement. 
Comparing with the actual time scale we see that the pho- 
toionisation time scale is relatively more important close to 
the source than further away due to shielding effect of the 
density enhancement. As time increases the shielding effect 
decreases and recombination time scale becomes more uni- 
form in the spatial domain. We note that these time scales 
can vary over several orders of magnitude across the ion- 
ising region so that it is not possible to describe the time 
dependent evolution using average recombination and pho- 
toionisation time scales. 

3.3 The role of time dependent term 

The calculations that we have presented so far have ne- 
glected the time dependent term ^^^^j-^ in the transfer 
equation. To illustrate the role played by this term, we have 
repeated the calculations for our model Model bS discussed 
in the previous section including this term. The results are 
presented in Figure ITUl 

The term with the time derivative in the transfer equa- 
tion 13 may be expected to play an important role when 
X >> ct. The neglect of this term could result in faster 
than light propagation of ionising fronts. The light travel 
time across our computational region is 6.5 x 10~^ Ma, and 
in our calculations, this condition is satisfied only at very 
early times which we have not attempted to resolve. The 
results of Figure [101 however, show that when taken in con- 
junction with the rate equations, this term has an additional 
effect on the state populations at all spatial locations yield- 
ing values that are lower than when the term is neglected. 
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Figure 7. The rate of change of the neutral density | \ on a logarithmic scale for the models with density 
enhancements; Model a (top left), Model b (top right), Model c (bottom left), and Model d (bottom right) 
showing the propagation of the ionisation fronts through a Gaussian density enhancement. The speed of the 
front is reduced and a time delay introduced due to the encounter. 




Figure 8. The spatial and time dependence of the recombination time scale (left) and the total time scale 
No/\i§^\ (right) for Model b3 which has a Gaussian density enhancement at a; = lOkpc 
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Figure 9. Tlie spatial and time evolution of the logarithm of 
the intensity of radiation for the models in Figure |H] Apart for a 
scaling factor, the inverse of the intensity is proportional to the 
photoionisation time scale 



The effect arises from tlie non-linearity of the rate equa- 
tions, and its magnitude depends on spatial location and 
time. Although we have neglected this term in our present 
models, our method allows for its inclusion with little addi- 
tional cost in computing time. This term may play a role in 
cosmological calculations. 

The approach to the steady state solution of this par- 
ticular model was shown in Figure |5| It is clear that none 
of the models in Figure |S| have reached a steady state in 
tmax 4 Ma, the assumed lifetime of the source indicating 
that full time dependent calculations are required to esti- 
mate the scale of ionisation domains. 




Figure 10. The spatial and time evolution of the logarithm of 
the difference between the calculated neutral density with and 
without the term ^dl/dt in the transfer equation for Model b3 
(see Figure|SJ. The solid line refers to .03 Ma, the broken curve to 
2 Ma and the dotted curve to 4 Ma. The light travel time across 
the computational domain is 6.5 X 10^'^ Ma. 



3.4 The global evolution of ionisation with time 

A good global indicator of the time evolution of ionisation 
can be obtained by plotting the number of neutral hydro- 
gen atoms in our fixed computational volume as a function 
of time t. This is shown in Figure [TTI for the three uniform 
density models that have been discussed. The time evolu- 
tion of the optical depths of these configurations can also be 
obtained from these diagrams by simply shifting the vertical 
scale by log ^ ~ 18. 

In the very low density model (Figure ITTji . the column 
density of neutral gas within the computational region is 
essentially constant decreasing only very slowly with time. 
For this model, the simple estimates based on eauations il5l 
and 1161 give tree ~ 630 Ma, fphot ~ 0.32 Ma. As the density 
is increased by a factor 10 and optical depth effects become 
more important, the column density of neutral gas is seen 
to decline more rapidly initially due to photoionisation, but 
then more slowly as recombination begins to play a signifi- 
cant role. For this model t^ac ~ 63 Ma. These effects are seen 
even more strongly in the highest density model considered 
for which t^cc '~ 0.63 Ma. 

We have not shown the corresponding curves for the 
models with density enhancements discussed in Figures |S| 
since they are identical to the curve with A^amb = 10~* cm~^ 
in Figure [TTI This is because all these four models have the 
same total column density in the computational region. The 
main effect of the introduction of the density enhancement 
is to change the spatial distribution of ionised gas, and the 
time evolution of the ionisation fronts, but not the total 
column density of ionised gas. 

The column density of ionised gas, and its evolution 
with time is therefore a strong function of the total column 
density and the incident radiation field and is not sensitive 
to the distribution of mass within the computational region. 

A related quantity that can be derived from our calcu- 
lations is the time evolution of the ratio of the total number 
of ionising photons in our computational region to the total 
number of ionised particles. We have derived this quantity 
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Figure 11. The time evolution of the logarithm of the neutral gas 
density within the entire computational volume for the uniform 
density Models a (bottom), b (middle) and c (top). This figure 
shows the effects of increasing the value of by a factor 100. 

^phot 

for our model with a central density enhancement (Model 
b4) and the results are shown in Figures [T^ 

Initially, the gas is neutral, and we therefore expect 9 to 
be large. As time proceeds, photons are removed from this 
region generating ionised particles. As photoionisations pro- 
ceeed, 6 declines rapidly, reaches a plateau, and rises again 
as the competing effect of recombination comes into play. 

3.5 Effects of increasing the number of continuum 
states 

The models that we have presented so far have been based 
on a model atom with one bound state and one continuum 
state. To allow for different spectral slopes for the source of 
ionising photons and a better description of the photoioni- 
sation cross section we need to include more than one con- 
tinuum state. 

We show in Figure [T^ the effect of including three con- 
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Figure 12. The time evolution of the logarithm of the ratio 8{t) 
of the total number of ionising photons to the total number of 
ionised particles in our computational box for Model b4. 



tinuum levels with A = 911, 422 and 196 A in Model b dis- 
cussed in Figure 131 b. The incident spectrum is assumed to 
be flat in wavelength but have the same total average in- 
tensity lirr ~ 10^ erg cm~^ s~^. The spatial profiles of 
the neutral density component and its time evolution are 
markedly different even though the maximum degree of ion- 
isation remains the same. We see that the central regions 
ionise more quickly and that a well defined front that prop- 
agates outwards can be seen within our computational re- 
gion. The intensity of radiation (Figures 1131 b, c and d) is 
now strongly frequency dependent reflecting the depen- 
dence of the photoionisation cross-section, with the high- 
est frequency showing the strongest attenuation close to the 
source at early times. A distant observer will thus see the 
source first at high frequencies, and then at lower frequencies 
as the Lyman limit is reached. 



3.6 Sensitivity to input parameters 

The ionisation state of the intergalactic medium depends on 
the nature of the sources of ionising radiation (modelled by 
lirr), their distribution in space and time, and the density 
structure of the ambient intergalactic medium. Our calcula- 
tions show that small uncertainties in these input parame- 
ters can have large effects on our conclusions on the state of 
ionisation of the intergalactic medium. 

To illustrate the sensitivity of our results on the state 
of ionisation to input parameters, we have repeated the cal- 
culations for Model b4, first by keeping lirr fixed, and de- 
creasing the ambient density by 10 percent, and then by 
keeping the ambient density fixed, and decreasing lirr by 
10 percent. The results of these two sets of calculations are 
shown in Figure 1141 A 10 percent variation in either the 
density, or the incident intensity translates to a a factor 3 in 
the neutral density on a time scale of 4 million years. These 
results show that the details of how re-ionisation proceeds 
requires accurate calculation of the density structure in the 
early universe, such as those that can be obtained from hy- 
drodynamic simulations (Abel et al. (1999)). 
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Figure 13. (a) (top left): The time evolution of the logarithm of 
the number density of neutral hydrogen at three different times 
t = 0,2,4 Ma. The parameters are as for the standard Model 
b except that three continuum levels (A = 911,422,196 A) are 
included: (b, c, d) (top right, bottom left and right): The spatial 
distribution and time evolution of the intensity of radiation at 
the the three wavelengths A = 911, 422, 196 A) respectively 



-4.5 



tji -5 . 5; 




12 3 4 

t [Ma] 




12 3 4 

t [Ma] 



Figure 14. The sensitivity of the neutral density to a reduction 
in the ambient density by 10% (upper panel) and to a reduction 
in the incident intensity by 10% (lower panel) for Model b4. 



4 MODEL RESULTS FOR MULTIPLE 
SOURCES 

In this section we present results of one set of calculations 
that were carried out with two sources using the generalisa- 
tion of the numerical method described in section 2.2 (Model 
si in Table 1). In these models, the first source is located 
a.t X — 0, turns on at t = 0, and shines for 5 Ma after which 
time it shuts ofT. A second source located at a; = 20 kpc 
turns on at a subsequent time t = 15 Ma, shines for an- 
other 5 Ma, and then shuts off. Each source has lirr ~ 10^ 
erg cm~^ s~^ and we assume a density distribution with a 
density enhancement centered at a;s = as in Model b4. 

We show in Figure IT^ the time evolution of the loga- 
rithm of the number density of neutral hydrogen. Initially, 
the behaviour is similar to the single source case, with the 
ionisation front propagating outwards with time starting at 
X = with the number density of neutral hydrogen decreas- 
ing with time at any given spatial position. Once the first 
source of radiation turns off (f = 5 Ma) , the gas begins to re- 
combine slowly, and the number density of neutral hydrogen 
increases at any given spatial location, but not significantly 
in the intervening 15 Ma. The second source then turns on 
and ionisation proceeds again, but starting from a medium 
which is already partially ionised. In the next 5 Ma the de- 
gree of ionsiation drops rapidly reaching values as low as 
~ 10~^ over most of the region at the end of the calcula- 
tion. 

The evolution of mean intensity J in our computational 
volume is shown in 115b . If we were located at a; = 10 kpc, 
mid way between the two sources, we would first see the 
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Figure 15. (a) The time evolution of the logarithm of the num- 
ber density of neutral hydrogen at five different times t = 0.6 Ma 
(solid), t = 5 Ma (long dashed), t = 10 Ma (short dashed), t = 
15 Ma (dashed - dotted) and t = 20 Ma (dotted) for the two 
source model (Model si), (b) The mean intensity in (x,t) space 
for Model si. 

radiation from the first source rise to reach a peak at t = 
5 Ma, and the source will then vanish from view almost 
instantaneously on a time scale corresponding to the light 
travel time from the source. The second source will then 
come into view on a similar time scale, and brighten slowly 
as the medium becomes significantly ionised again. 

The net effect of having two sources of ionisation that 
turn on at different times is to ionise a significantly larger 
volume of gas to a higher degree of ionisation. These calcu- 
lations demonstrate that after the first sources turn on, the 
universe will ionise very quickly as subsequent sources turn 
on, even though the lifetime of the sources of radiation may 
be much shorter than the time between the sources. 



5 DISCUSSION AND CONCLUSIONS 

We have presented a robust numerical scheme that is un- 
conditionally stable and allows the calculation of the time 
dependent evolution of ionisation domains with full time and 
spatial resolution. The method allows for the explicit inclu- 
sion of the time dependence of the radiative intensity in the 
transfer equation. 

Our illustrative 1-D examples were chosen to mimic con- 
ditions in the early universe and highlight the basic processes 
that are involved in the development of ionisation domains 
around UV sources of radiation. The physics of the non- 
LTE ionisation process is of course well understood, and 
has been incorporated at different levels of approximation 



in many previous investigations of cosmological reionisation. 
As a source first turns on, and illuminates a neutral medium, 
photoionisations will dominate everywhere except that they 
will be less effective in the outer regions because of the larger 
optical depth towards the illuminating source. The neutral 
density will therefore initially decline at different rates at 
different spatial locations in the medium with the neutral 
density in the inner regions declining more rapidly. As the 
neutral density decreases, recombinations become relatively 
more important and slows down the rate of decline. In addi- 
tion the ionisation of the regions closest to the source makes 
this gas more transparent enhancing the photoionisation of 
adjacent outer regions. These two processes combine to yield 
an expanding ionisation 'front' that separates the mainly 
ionised and mainly neutral gas. The front expands at a rate 
that is determined by the local density and the intensity Urr 
of the source of irradiation. As t increases, the amplitude of 
^^f- 1 decreases, but an equilibrium situation is not reached 
in any of our calculations within the lifetime of the source 
(~ 4 Ma). 

The above simple picture can change significantly when 
density fiuctuations are taken into consideration. We have 
investigated the effects that come into play when an outward 
propagating ionising front encounters a Gaussian density en- 
hancement of given peak density and half width. The density 
enhancement slows down the speed of the front in a com- 
plex manner depending on both of these parameters. Before 
the front encounters the enhancement, it has a speed that 
is nearly constant and equal to the value appropriate to the 
ambient mean density and the input irradiating intensity 
lirr- As the front emerges from the density enhancement, 
the speed reaches a value that is generally lower than be- 
fore the encounter. The density enhancement 'delays' the 
propagation of the ionisation front by a time At/ that is 
strongly dependent on the peak density and the half width 
of the assumed Gaussain density enhancement. Our calcu- 
lations have shown that delays of the order of a few tens of 
million years are possible even for modest density enhance- 
ments. The delay could therefore exceed the lifetime of the 
ionising source, so that the density enhancement in effect 
stops the propagation of the front. Due to the long recombi- 
nation times, such regions will remain ionised, enshrouded 
by neutral gas, until the next source of radiation turns on 
in its vicinity and continues the process of ionsiation. Our 
calculations have shown that even though the lifetime of a 
source of ionising radiation may be much smaller than the 
time between the emergence of such sources, ionisation will 
proceed very effectively once started. 

Various generalisations, such as extension to three di- 
mensions are easily implemented. Helium ionisation can be 
readily incorporated by the use of more levels (and there- 
fore more frequency points) and the cosmological expansion 
term by implicitly discretizing the frequency variable. Our 
future plans are to apply our numerical method to outputs 
of cosmological simulations. 
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